
\chapter{蒙特卡洛方法}

我个人相信，人类的博彩业可以追溯到史前文明，每当色(shai三声)子被掷下，金钱就从一个人流向另一个人。
在欧洲的摩纳哥公国，有一个著名的城市名叫蒙特卡洛，每天都在进行大量的赌博活动。
无数个色子在光洁的桌面上滚动，伴随着男男女女的欢呼与悲叹。
光洁的色子上，镌刻着人类命运的点数，随机地倒向阴暗与光明。
随机性，是那么让人类痴迷。
于是聪明绝顶的冯洛伊曼，决定用蒙特卡洛这个神秘的名字，命名一类用到随机数的方法。

\subsection{随机数的期望值与方差}
如果一个随机变量 $X$ 的概率密度为 $f(x)$（不失一般性，离散分布的随机数可以用带有 $\delta$ 函数的 $f(x)$ 描述），则其期望值为
\begin{eqnarray}
E(X) = \int^\infty_{-\infty} x f(x) dx,
\end{eqnarray}
其方差为
\begin{eqnarray}
D(X) = \int^\infty_{-\infty} ( x - E(X) )^2 f(x) dx.
\end{eqnarray}

如果相互独立的随机变量 $X_1, X_2, \cdots, X_n$ 的概率密度分别为 $f_1(x), f_2(x), \cdots, f_n(x)$，那么 $X_1 + \cdots + X_n$ 的期望和方差分别是多少？
\begin{eqnarray}
E(X_1 + \cdots X_n) &=& \int dx_1 \int dx_2 \cdots \int dx_n f_1(x_1) f_2(x_2) \cdots f_n(x_n) (x_1 + x_2 +\cdots x_n) \nonumber\\
&=& \sum^n_{i=1} \int dx_1 \int dx_2 \cdots \int dx_n f_1(x_1) f_2(x_2) \cdots f_n(x_n) x_i \nonumber\\
&=& E(X_1) + \cdots E(X_n) \nonumber\\
D(X_1 + \cdots X_n) &=& \int dx_1 \int dx_2 \cdots \int dx_n (x_1 + x_2 + \cdots + x_n - E(X_1 + \cdots + X_n))^2 \nonumber\\
															&& \times f_1(x_1) f_2(x_2) \cdots f_n(x_n) dx_1 dx_2 \cdots dx_n \nonumber\\
												&=& \int dx_1 \int dx_2 \cdots \int dx_n \left[ (x_1-E(X_1)) + (x_2 - E(X_2)) + \cdots + (x_n - X_n) \right]^2 \nonumber\\
															&& \times f_1(x_1) f_2(x_2) \cdots f_n(x_n) dx_1 dx_2 \cdots dx_n \nonumber\\
												&=& D(X_1) + D(X_2) + \cdots D(X_n)
\label{EDsum}
\end{eqnarray}

如果一个随机变量 $X$ 的概率密度函数为 $f(x)$，$\lambda$ 为一个正实数，那么 $\lambda X$ 的概率密度、期望与方差分别是多少？
\begin{eqnarray}
P\left\{ \lambda X \leq \alpha \right\} = P\left\{ X \leq \frac{\alpha}{\lambda} \right\} = \int^\frac{\alpha}{\lambda}_{-\infty} f(x) dx,
\end{eqnarray}
所以 $\lambda X$ 的概率密度函数为
\begin{eqnarray}
g(\alpha) = \frac{d}{d\alpha} P\left\{ \lambda X \leq \alpha \right\} = \frac{d}{d\alpha} \int^\frac{\alpha}{\lambda}_{-\infty} f(x) dx = \frac{1}{\lambda} f(\alpha/\lambda).
\end{eqnarray}
$\lambda X$ 的期望值为
\begin{eqnarray}
E(\lambda X) = \int^\infty_{-\infty} \frac{\alpha}{\lambda} f(\alpha/\lambda) d\alpha = \lambda E(X),
\label{ElambdaX}
\end{eqnarray}
$\lambda X$ 的方差为
\begin{eqnarray}
D(\lambda X) &=& \int^\infty_{-\infty} \left[\alpha - E(\lambda X)\right]^2 \frac{1}{\lambda} f(\alpha/\lambda) d\alpha \nonumber\\
				&=& \int^\infty_{-\infty} \left[\alpha - \lambda E(X) \right]^2 \frac{1}{\lambda} f(\alpha/\lambda) d\alpha \nonumber\\
				&=& \lambda^2 \int^\infty_{-\infty} \left[ \frac{\alpha}{\lambda} - E(X)\right]^2 f(\alpha/\lambda) d \frac{\alpha}{\lambda} \nonumber\\
				&=& \lambda^2 D(X)
\label{DlambdaX}
\end{eqnarray}

\subsection{大数定理与中心极限定理}
蒙特卡洛方法的理论依据是概率论中的大数定理与中心极限定理。
根据大数定理，随机抽样的统计平均会逼近理论预期；而中心极限定理可以用来估算抽样结果的误差。
所以这里简述这两个定理的证明，如果需要更多细节，请参考\cite{ProbabilityTheory}。

\subsubsection{切比雪夫（Chebyshev）不等式}

设随机变量 $X$ 的期望 $E(x)$ 和 方差 $D(x)$ 都存在，则 $\forall \epsilon>0$，都有切比雪夫不等式：
\begin{eqnarray}
P\left\{ | X - E(X) | \leq \epsilon \right\} \geq 1- \frac{ D(X) }{ \epsilon^2}
\end{eqnarray}
证明：假设概率密度函数为 $f(x)$，则有
\begin{eqnarray}
P\left\{ | X - E(X) | \geq \epsilon \right\} &=& \int_{ | X - E(X) | \geq \epsilon} f(x) dx \nonumber\\
				&\leq& \int_{ | X - E(X) | \geq \epsilon} \frac{ \left[x- E(X)\right]^2 }{\epsilon^2} f(x) dx \nonumber\\
				&\leq& \frac{1}{\epsilon^2} \int^\infty_{-\infty} (x- E(X))^2 f(x) dx \nonumber\\
				&=& \frac{ D(X) }{ \epsilon^2}.
\end{eqnarray}
所以
\begin{eqnarray}
P\left\{ | X - E(X) | \leq \epsilon \right\} = 1- P\left\{ | X - E(X) | \geq \epsilon \right\}
 \geq 1- \frac{ D(X) }{ \epsilon^2}，
\end{eqnarray}
即证明了切比雪夫不等式。

\subsubsection{切比雪夫大数定理}

设独立随机变量 $X_1, X_2, \cdots, X_n, \cdots$ 具有相同的数学期望和方差，$E(X_i)=\mu, D(X_i) = \sigma^2 (i=1,2,\cdots)$，则 $\forall \epsilon$，有
\begin{eqnarray}
\lim_{n \rightarrow \infty} P \left\{ | \frac{1}{n} \sum\limits^n_{i=1} X_i - \mu | \leq \epsilon \right\} =1.
\label{large_number_theorem}
\end{eqnarray}
证明：根据公式（\ref{EDsum}-\ref{DlambdaX}），
\begin{eqnarray}
E( \frac{1}{n} \sum \limits^n_{i=1} X_i ) &=& \mu, \nonumber\\
D( \frac{1}{n} \sum \limits^n_{i=1} X_i) &=& \frac{ \sigma^2}{n}
\end{eqnarray}
所以，根据切比雪夫不等式，$\forall \epsilon >0$,
\begin{eqnarray}
P\left\{ | \frac{1}{n} \sum\limits^n_{i=1} X_i - \mu | \leq \epsilon \right\} \geq 1- \frac{ \sigma^2 }{ n \epsilon^2},
\end{eqnarray}
所以
\begin{eqnarray}
\lim_{n \rightarrow \infty} P \left\{ | \frac{1}{n} \sum\limits^n_{i=1} X_i - \mu | \leq \epsilon \right\}
\geq \lim_{n \rightarrow \infty} \left[1- \frac{ \sigma^2 }{ n \epsilon^2} \right]
 =1.
\end{eqnarray}
这个数不可能超过1，所以它的极限必等于1。这样，我们就证明了切比雪夫大数定理。

在实践中，我们每次抽样的结果都是随机的，有一个概率分布和确定的期望、方差；切比雪夫大数定理意味着，随着抽样次数增多，抽样结果的统计平均必然收敛于每次抽样的期望值。

\subsubsection{中心极限定理}
{\bf 独立同分布的中心极限定理}：假设随机变量 $X_1, \cdots X_n$ 相互独立且遵循同一分布，期望值为 $\mu$，方差为 $\sigma^2$，则
\begin{eqnarray}
\lim_{n \rightarrow \infty} \frac{ \sum\limits^n_{i=1} X_i - n \mu }{ \sqrt{n} \sigma } \sim N(0,1),
\label{central_limit_theorem}
\end{eqnarray}
其中 $N(0,1)$ 表示期望为 $0$，方差为 $1$ 的正态分布。这是一个美妙的定理，因为它不依赖于 $X_i$ 的具体概率密度函数。它说明正态分布是大自然的选择。

为了证明这个定理，我们先简单介绍概率论中所谓“特征函数”：如果一个随机变量 $X$ 的概率密度函数为 $f(x)$，则其特征函数为 $\varphi_X (t) = \int^\infty_{-\infty} e^{i t x } f(x) dx$，即正比于概率密度函数的傅里叶变换。
如果一个随机变量的概率密度函数确定了，其特征函数也唯一确定，反之亦然，所以特征函数也反映了随机变量的性质。
例如，特征函数 $\varphi_X(t)$ 的一阶导数、二阶导数在 $t=0$ 处的值反映了随机变量 $X$ 的期望和方差。
\begin{eqnarray}
\frac{d}{dt} \varphi_X(t)|_{t=0} &=& \frac{d}{dt} \int^\infty_{-\infty} e^{i t x} f(x) dx |_{t=0} \nonumber\\
																	&=& \int^\infty_{-\infty} i x f(x) dx = i E(X) \nonumber\\
\frac{d^2}{dt^2} \varphi_X(t)|_{t=0} &=& - \int^\infty_{-\infty} x^2 f(x) dx = - D(X) - E(X)^2
\label{characteristic_derivative}
\end{eqnarray} 
正态分布 $N(0,1)$ 的概率密度函数为
\begin{eqnarray}
f_{N(0,1)}(x) = \frac{1}{ \sqrt{2\pi} } e^{-\frac{x^2}{2}},
\end{eqnarray}
其特征函数为
\begin{eqnarray}
\varphi_{N(0,1)}(t) &=& \int^\infty_{-\infty} \frac{1}{\sqrt{2\pi}} e^{ - \frac{x^2}{2} + itx}dx = e^{-\frac{1}{2} t^2}.
\label{characteristic_normal_distribution}
\end{eqnarray}
最后一步用了泊松积分，如果你感兴趣，可以参考\cite{Calculus4}第123页。
值得注意，$f_{N(0,1)}(x)$ 与 $\varphi_{N(0,1)}(t)$ 具有相同的形式，具有对称的美。

此外，如果 $X_1, X_2$ 是相互独立的随机变量，其概率密度函数分别为 $f(x), g(x)$，则 $X_1 + X_2$ 的概率密度函数为
\begin{eqnarray}
h(x) &=& \frac{d}{dx} \int^\infty_{-\infty} d \xi f(\xi)  \int^{ x- \xi}_{-\infty} g(\eta) d\eta = \int^\infty_{-\infty} f(\xi) g(x-\xi) d\xi,
\end{eqnarray}
其傅里叶展开为
\begin{eqnarray}
\varphi_{X_1 + X_2} (t) = \int^\infty_{-\infty} \int^\infty_{-\infty} f(\xi) g(x-\xi) d\xi e^{itx} dx = \varphi_{X_1}(t) \varphi_{X_2}(t).
\end{eqnarray}
这个结论很容易推广到 $n$ 个独立随机变量的和，即
\begin{eqnarray}
\varphi_{X_1 + \cdots X_n} = \varphi_{X_1}(t) \cdots \varphi_{X_n}(t).
\label{characteristic_sum}
\end{eqnarray}

下面我们证明独立同分布中心极限定理，即公式（\ref{central_limit_theorem}），这个思路来自\cite{wiki}中``central limit theorem”词条，略有修改。
记
\begin{eqnarray}
Z_n = \frac{ \sum\limits^n_{i=1} X_i - n \mu }{ \sqrt{n} \sigma } = \sum\limits^n_{i=1} Y_i,
\end{eqnarray}
其中 $Y_i \equiv  \frac{1}{\sqrt{n}} \frac{X_i - \mu}{\sigma}$，根据公式（\ref{EDsum}-\ref{DlambdaX}），显然 $E(Y_i)=0, D(Y_i)=\frac{1}{n}$。
根据公式（\ref{characteristic_sum}），$Z_n$ 的特征函数为
\begin{eqnarray}
\varphi_{Z_n}(t) = \varphi_{Y_1}(t) \cdots \varphi_{Y_n}(t) = \varphi_{Y_1}(t) ^n.
\end{eqnarray}
根据公式（\ref{characteristic_derivative}），我们有
\begin{eqnarray}
\varphi_{Z_n}(t) = \varphi_{Y_1}(t)^n &=& \left\{ 1 + iE(Y_1)t - \frac{1}{2} \left[ D(Y_1) + E(Y_1)^2 \right] t^2 + \cdots \right\}^n \nonumber\\
									&=& \left\{ 1 - \frac{1}{2n} t^2 + \cdots \right\}^n
\end{eqnarray}
故 $n\rightarrow \infty$ 时，有
\begin{eqnarray}
\lim_{n \rightarrow \infty} \varphi_{Z_n}(t) &=& \lim_{n \rightarrow \infty}  \left\{ 1 - \frac{1}{2n} t^2 + \cdots \right\}^n \nonumber\\
																							&=& e^{-\frac{1}{2}t^2}
\end{eqnarray}
比较上式与公式（\ref{characteristic_normal_distribution}），可知
\begin{eqnarray}
\lim_{n \rightarrow \infty} \frac{ \sum\limits^n_{i=1} X_i - n \mu }{ \sqrt{n} \sigma } \sim N(0,1),
\end{eqnarray}
这样我们便证明了独立同分布的中心极限定理。

\subsection{蒙特卡洛方法}
\label{Monte-Carlo}
\subsubsection{蒙特卡洛方法的原理}
在蒙特卡洛方法中，如果我们需求物理量 $\mu$，可以构造算法，使得随机变量 $X$ 的期望值 $E(X)=\mu$，然后进行随机数值试验，对 $X$ 进行随机抽样。
根据大数定理，即公式（\ref{large_number_theorem}），有
\begin{eqnarray}
\lim_{n \rightarrow \infty} P \left\{ \frac{1}{n} \sum\limits^n_{i=1} X_i \right\} = \mu,
\end{eqnarray}
说明在抽样次数趋于无穷时，抽样平均一定收敛于所求物理量 $\mu$，所以抽样平均可以作为 $\mu$ 的数值近似。

在实践中我们的随机抽样次数永远是有限的，抽样结果的统计平均并没有精确地收敛于所求物理量，所以需要估算误差。
而中心极限定理（公式（\ref{central_limit_theorem}））告诉我们，
\begin{eqnarray}
\lim_{n \rightarrow \infty} \frac{ \frac{1}{n} \sum \limits^n_{i=1} X_i - \mu }{ \sigma/ \sqrt{n} } \sim N(0,1),
\end{eqnarray}
它说明，进行大样本抽样时，抽样结果的统计平均近似地服从高斯分布，所以我们可以根据高斯分布，定出与{\bf 置信度}相对应的 {\bf 置信区间}，即数值误差范围。
对于任意给定的置信度 $1-\alpha$（$0<\alpha<1$），有
\begin{eqnarray}
\lim_{n \rightarrow \infty} P \left\{ |\frac{1}{n} \sum \limits^n_{i=1} X_i - \mu| < \frac{X_\alpha \sigma}{\sqrt{n}} \right\} = \frac{2}{\sqrt{2\pi}} \int^{X_\alpha}_0 e^{-x^2/2} dx = 1-\alpha,
\end{eqnarray}
即 $(\frac{1}{n} \sum \limits^n_{i=1} X_i - \frac{X_\alpha \sigma}{\sqrt{n}}, \frac{1}{n} \sum \limits^n_{i=1} X_i + \frac{X_\alpha \sigma}{\sqrt{n}})$ 即为与置信度 $1-\alpha$ 相对应的置信区间。
我们得到的结论是：物理量 $\mu$ 的值有 $1-\alpha$ 的可能性存在于这个区间内。
$\sigma$ 可以由
\begin{eqnarray}
\sigma \approx \sqrt{\frac{1}{n} \sum \limits^n_{i=1} X^2_i - (\frac{1}{n} \sum \limits^n_{i=1} X_i)^2 }
\end{eqnarray}
来近似地估计。

\subsubsection{“随机数”的产生}
量子力学中的测量在哥本哈根阐释下具有随机性（几率阐释），但是计算机不可能让人类猜大小，它们遵从人类的指令，每步给出确定的结果。
要用计算机模拟随机数，只好设法生成一个毫无规则的数列，其中的数在定义域中任意处取值的机会均等；每次要在计算中“生成”一个随机数，便从这个无规则的数列中取用一个。
数值计算和模拟中用的“随机数”都是这样的无规数列，我们称之为“伪随机数”。
产生这种数列的算法确定以后，给它一个初始数字（即“种子”），它便可以生成任意多的随机数；如果种子确定，后面的每个数便都确定了。
为了在每次运行程序时都产生不同的无规数列，往往将“种子”设作计算时计算机上的时间（从1970年到计算的那一刻的总秒数），这样，每次运行程序时随机数种子都不同，所以“伪随机数”更具有一般性。

伪随机数要满足一定的要求\cite{KZYan}：1. 任何算法下，伪随机数都有周期，经过一个周期以后会出现重复的数字，所以伪随机数的周期必须足够长；
2. 有良好的随机性，否则不能模拟随机过程，一个检验方法是在考察 $x_n$ 与 $x_{n+i}$ 的关联，如果出现关联性，就达不到要求；
3. 计算伪随机数的速度必须快，否则不能进行大量抽样。

此处只介绍“乘加同余法”，
\begin{eqnarray}
I_{j+1} = a I_j + c (\rm mod~~ m),
\end{eqnarray}
其中 $a,c,m$ 均为正整数，产生概率密度为常数的随机整数，当 $m$ 为很大的整数时具有较大的周期。
其他概率密度函数的随机数可以由此出发构造出来。

\subsubsection{多重数值积分}
前文中提到过，多重积分的计算量与积分重数是指数关系，而且算法精度会下降，所以对于重数较高的数值积分，很难用传统的方法计算。
根据\ref{Monte-Carlo}中的分析，我们可以用蒙特卡洛方法进行数值估计。
我们可以将被积函数在积分区间内的平均值作为所求物理量，在区域内随机取点，计算函数值并做统计平均，得到
\begin{eqnarray}
\int_V f dV &\approx& \langle f \rangle V \pm V X_\alpha \frac{\sigma}{ \sqrt{n} } \nonumber\\
						&\approx& \langle f \rangle V \pm V X_\alpha \frac{ \sqrt{ \frac{1}{n} \sum \limits^n_{i=1} f^2(\vec{x}_i) -  \left[ \frac{1}{n} \sum \limits^n_{i=1} f(\vec{x}_i) \right]^2 } }{ \sqrt{n} },
\label{multi-dimension-integral}
\end{eqnarray}
作为置信度为 $1-\alpha$ 的数值估计，其中 $\frac{2}{\sqrt{2\pi}} \int^{X_\alpha}_0 e^{-x^2/2} dx = 1 - \alpha$。

公式（\ref{multi-dimension-integral}）中最后一项分数线上面是 $\sigma$ 的估计，随着 $n$ 的增大理应收敛于一个稳定的常数。
所以蒙特卡洛多重积分的数值误差 $\sim (\frac{1}{n})^{1/2} $，前文中提到，如果一个普通数值积分算法（如均值法、梯形法、辛普生方法）的误差在一维情况下 $\sim (\frac{1}{n})^{a}$，则在 $d$ 维情况下误差为 $\sim (\frac{1}{n})^{a/d}$。所以，当 $d$ 很大时，有可能 $a/d << 1/2$，蒙特卡洛方法的威力会体现出来。

对于振荡剧烈的被积函数，如果在积分区域内均匀抽样，$\sqrt{ \frac{1}{n} \sum \limits^n_{i=1} f^2(\vec{x}_i) -  \left[ \frac{1}{n} \sum \limits^n_{i=1} f(\vec{x}_i) \right]^2 }$ 的值可能很大，很难给出置信度大，而置信区间小的结果。
可以按照概率密度函数 $g(\vec{x})$ 取抽样点，即 $g(\vec{x})$ 的值较大的地方，取点密度更大，$g(\vec{x})$ 的值较小的地方，取点密度更小。
相应地，对所有抽样结果乘以 $\frac{1}{g(\vec{x})}$ 的权重，再做平均。
可以证明，这样做不影响对 $f(x)$ 或 $f^2(x)$ 的估计，却有可能显著缩小抽样值的方差，在 $1-\alpha$ 的置信度下相应置信区间为
\begin{eqnarray}
\int_V f dV \approx V \langle f/g \rangle \pm V X_\alpha \sqrt{ \frac{ \langle f^2/g \rangle - \langle f/g \rangle^2 }{ n }},
\label{multi-dimension-integral2}
\end{eqnarray}
其中 $n$ 为抽样次数。
所以选取适当的函数 $g(\vec{x})$，即可巧妙地缩小蒙特卡洛估计的置信区间。
可以证明，当 $g(\vec{x}) = \frac{f(\vec{x})}{ \langle f/g \rangle }$ 时，公式（\ref{multi-dimension-integral2}）中的置信区间为零。

\subsubsection{麦氏模拟（Metropolis）}
对于正则系综力学量在不同微观状态中的统计平均值为
\begin{eqnarray}
\langle A \rangle = \frac{ \int A e^{-\beta E} d \Omega }{ \int e^{-\beta E} d \Omega}
\end{eqnarray}
如果不同微观状态的抽样概率正比于 $e^{-\beta E}$，其中 $E$ 为微观状态对应的体系总能量，则
\begin{eqnarray}
\langle A \rangle = \frac{1}{n} \sum\limits^{n}_{i=1} A_i
\end{eqnarray}

{\bf 马尔科夫（Markov）链}

如果能构造一个链，链上每一点表示一个微观状态，丢掉链的前面若干个点后，剩余的点满足正则分布——即任意一个微观状态 $s$ 出现的概率正比于 $e^{-\beta E_s}$——则从链上抽样，即可实现上述麦氏模拟。

怎样构造马尔科夫链呢？链上每两个相邻的微观状态之间由一个微小变化连接。如果链上每一点从状态 $r$ 变为状态 $s$ 的概率为 $\omega(\vec{x}_r - \vec{x}_s)$，反之为 $\omega(\vec{x}_s - \vec{x}_r)$，则一段时间后
\begin{eqnarray}
\Delta N_{r \rightarrow s} &=& N_r \omega( \vec{x}_r - \vec{x}_s ) - N_s \omega(\vec{x}_s - \vec{x}_r) \nonumber\\
														&=& N_r \omega( \vec{x}_s - \vec{x}_r ) \left[ \frac{ \omega(\vec{x}_r - \vec{x}_s) }{ \omega(\vec{x}_s - \vec{x}_r)} - \frac{ N_s }{ N_r } \right]
\end{eqnarray}
如果算法确保
\begin{eqnarray}
\frac{ \omega( \vec{x}_r - \vec{x}_s ) }{ \omega( \vec{x}_s - \vec{x}_r) } = \frac{ e^{-\beta E_s} }{ e^{-\beta E_r}},
\end{eqnarray}
则足够长的时间后，$N_s / N_r$ 必然收敛于 $e^{-\beta E_s} / e^{-\beta E_r}$。

所以只要 $\omega( \vec{x}_r - \vec{x}_s ) / \omega( \vec{x}_s - \vec{x}_r)  = e^{- \beta (E_s - E_r)}$，即可构造马尔科夫链，进行抽样，实现麦氏模拟。
例如，可以选取
\begin{eqnarray}
\omega( \vec{x}_r - \vec{x}_s ) = \left\{
\begin{aligned}
&exp \left\{ - \beta (E_s - E_r) \right\}, &E_s > E_r;\\
&1, ~~~~~~~~~~~~~~~~~~~~~~~~&E_s < E_r.
\end{aligned}
\right.
\label{omega}
\end{eqnarray}

\subsubsection{Ising模型的模拟}

二维Ising模型的模拟思路是这样的：首先定义$L\times L$的spinor点阵，即$n=L^2$个spinor呈正方形排布。
利用周期性边界条件，右边界上的spinor与左边界上的spinor相邻，而上边界的spinor与下边界上的spinor相邻，即所有spinor的地位都完全一样。
这样，用一个有限大的点阵模拟spinor晶体的整体性质。
在一定温度下，这个晶体中的所有spinor处在正则分布，即它们的取向、以及相关的相互作用可能处在各种微观状态，相应的概率满足正态分布。
所以，要模拟材料的性质，就需要对方阵中所有微观状态进行统计，得到系综平均。
而麦氏模拟可以实现这一点。

如上文中产生马尔科夫链，使得链上每个节点对应一个微观状态，不同微观状态之间的转变概率为(\ref{omega})式。
而这个转变定义为翻转方阵中任意一个spinor，比如第i个，相对应的前后微观状态能量之差为
\begin{equation}
E_s - E_r = 2 J \vec{S}_i \cdot \sum^{4}_{j=1} \vec{S}_j,
\label{DeltaE}
\end{equation}
其中，$\vec{S}_i$为第i个spinor翻转前的自旋，而$\vec{S}_j$是它的四个最近邻。
将(\ref{DeltaE})带入(\ref{omega})，即得到不同微观状态之间的转换概率。

所以，我们通过随机翻转方阵中的spinor，并控制翻转的概率，即可构造马尔科夫链。
舍弃最开始若干个微观状态，即可进行抽样，模拟正则分布下的spinor晶体。

%\section{Variational Monte Carlo}
